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ABSTRAGT 

We present iL-band observations of P Pic with the Gemini Planet Imager’s (GPFs) polarimetry 
mode that reveal the debris disk between ^ 0.3" (6 AU) and ^ 1.7" (33 AU), while simultaneously 
detecting P Pic h. The polarized disk image was fit with a dust density model combined with a Henyey- 
Greenstein scattering phase function. The best fit model indicates a disk inclined to the line of sight 
{p = 85.27°1o. 19) ^ position angle Opa = 30.35°lo.28 (slightly offset from the main outer disk, 

OpA ^ 29°), that extends from an inner disk radius of 23.6lo.6 to well outside GPFs field of view. 

In addition, we present an updated orbit for P Pic b based on new astrometric measurements taken in 
GPFs spectroscopic mode spanning 14 months. The planet has a semi-major axis of a = 9.2 ^q 4 AU, 
with an eccentricity e < 0.26. The position angle of the ascending node is Q = 31.75° ±0.15, offset 
from both the outer main disk and the inner disk seen in the GPI image. The orbital fit constrains 
the stellar mass of P Pic to 1.61 ±0.05Mq. Dynamical sculpting by p Pic b cannot easily account for 
the following three aspects of the inferred disk properties: 1 ) the modeled inner radius of the disk is 
farther out than expected if caused by P Pic 6 ; 2 ) the mutual inclination of the inner disk and P Pic b 
is ^ 4°, when it is expected to be closer to zero; and 3) the aspect ratio of the disk (ho = 0.137^0 ooe) 
is larger than expected from interactions with p Pic b or self-stirring by the disk’s parent bodies. 

Keywords: planet-disk interactions, techniques: polarimetric, astrometry, planets: individual {P Pic 
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1. INTRODUCTION 

The dynamical interactions between exoplanets and 
their local debris disks provide a unique window into the 
understanding of planetary system architectures and evo¬ 
lution. In this regard, the (3 Pic system is important as it 
is one of the rare cases where both a planet and a debris 
disk have been directly imaged. 


The 6 Pic s ystem first garnered interest after [Smith fc 
Terrile (1984) followed up a prominent Infrared A stro¬ 
nomical Satellite (IRAS) infrared excess detection (An- 


mann|19^ ) and imaged an edge-on circumstellar disk in 


dust scattered light. Since then, many observational and 
theoretical studies have helped to paint a picture of a dy- 


([Roberge et 


directly ima 

ged ^10-12 Mj planet ( 

Lagrange et al.||2009i 

Snellen et al. 

2014 Chilcote et al. 

I2U15D, an asymmet- 

ric debris disJ 

k (Lagage Pantin |] 

9941 IKalas & Jewitt 

19951), infalling small bodies (Beust iVlorbidelli 1996 

Kiefer et al. 

2014|), multiple planetesimal belts (jUkamoto 


aT 


and a circling gas cloud that 


indicate a recent collision between planetesimals (Dent 
et al.|[2M4 ). Here we examine the nature of the dynam- 


may 


ical relationship between the planet, [3 Pic 6, and the 
debris disk using polarimetric imaging and modeling of 
the innermost region of the disk. 

The overall structure of the disk—a depleted inner re¬ 
gion, an extended outer region, and an ap parent warp— 
has been w ell-established in the literature. ISmith fc Ter-I 
rile (1984) originally used optical depth arguments to 
inter that p Pic’s disk must be depleted of gra ins in¬ 


terior to a radius of ^ 30 AU. Burrows et al. ([1995 ) 
used HST/WFPC2 to image the disk in optical scat¬ 
tered light and described qualitatively a vertical warp 
in the midplane structure somewhere between 1.5" and 
10" radius. The first quantitative measurements of the 
midplane warp were derived from ground-based adap- 
tive optics (AO) obse rvations in the near infrared (NIR; 
Mouillet et al. 1997). In these data, the peak height 
of the warp is at ^ 3" radius, ^ 58 AU assumin g he¬ 
liocentric distance of 19.44 pc (van Leeuwen||2007|), and 
corresponds to 3° deviation frorn the position angle (PA) 
of the midplane measured beyond ^100 AU. 

Two geometrical interpretations of the apparent warp 
have been proposed. The first is that we are observing 
a single disk warped by forcing from a planet on an in¬ 
clined orbit . Using numerical mo dels and semi-analytic 
arguments, Mouillet et al. (1997) demonstrated that a 
planet inclined by 3''-5'' to a hypothetical disk can repli¬ 
cate the observed structure via a secular perturbation. 
The inferred mass of the planet depends on when the 
planet’s orbit was perturbed out of coplanarity, because 
in this paradigm the warp propagates radiahy outw ard 
on million year timescales. Augereau et al. (2001) ap¬ 
plied this model to explain several other observational 
features of the disk such as the larger scale asymmetries. 

Alternatively, the structure could be composed of two 
disks, with symmetric linear morphologies, superimposed 
on the sky plane. Two disks would appear to create a 
warp in the midplane of the primary disk because of a 
^ 3° difference in position angle. Based on high an¬ 
gular resolution optical data obtained wit h HST/STIS 
that clearly showed the warp component, [Heap et al. 


(2000) postulated that the sky plane contains “two disks 
o'" apart.” This interpretation is also favored in subse¬ 
quent studies based on multi-color HST/ACS/H RC ob¬ 
servations of p Pic’s disk (jGohmowski et al.|2006|). More 
detailed analytic modeling of these data are consistent 
with two dis ks with a relative po si tion angle on the sky 


of 3.2±1.3° (Ahmic et al. 2009). Ahmic et al. (2009) 


also find that the fainter inclined disk has a line 3f sight 
inclination 6.0±1.0°, whereas the brighter, primary disk 
is consisten t with being exactly edge-on. More recently 


Apai et al. (2015) presented a re-reduction of the early 
HbT'/bTlb observations, coupled with newer observa¬ 
tions obtained 15 years apart. They found that these 
observations were consistent with the two-disk interpre¬ 
tation, but they also examine a scenario where f3 Pic b is 
perturbing the disk. 

The perturbing planet scenario requires a planet with 
a mass, semi-major axis, and mutual inclination with re- 
spect to the fiat outer disk sufficient to create the warp. 


Lagrange et al. (2009) discovered [3 Pic 6, a planet with 
a mass and separation appropriate for creating the warp; 
with additional astrometric measurements, its or bit was 
constrained to a ^ 9 AU, i ^ 89° an d e ^ 0.1 (Chau- 


vin et al. 2012 Macintosh et al. 2014). If the planet is 
secularly perturbing the disk, we expect it to be in the 
same plane as the inner disk and misaligned from the flat 
outer disk (though it may appear to be aligned in projec¬ 
tion). One technical challenge is that the planet location, 
the inner warp and the outer disk have been measured 
on different angular scales and are detected using dif¬ 
ferent observing strategies. Therefore, systematic errors 
in the position angle calibrations between different data 
sets lead to uncertainty in the rel ative orientations of 
these three structures. For example, Currie et al. (2011) 
reported that t he planet’s orbit is mi saligned with the 
inner disk, but Lagrange et al. (2012) noted that they 
are consistent with alignment when all sources of error 
are accounted for. 


Lagrange et al. (2012) attempted to solve these prob- 
lems by constructing observations where a single instru¬ 
ment is used to simultaneously detect both the planet 
and the disk. The results show (3 Pic b positioned 2-4° 
above the southwest disk midplane at the 2010 epoch of 
observation (“above” means north of the SW midplane 
or at a larger PA than the SW midplane). Therefore, 
/3 Pic 6’s orbit is not coplanar with the main, flat, outer 
disk. Instead the position above the main midplane in 
the SW is in the direction of the warped component. 
This projected misalignment is consistent with the nec¬ 
essary mutual inclination between the planet’s orbit and 
the main, flat, outer disk. 

A different technical challenge is imaging the disk along 
the minor axis direction very close to the star in order 
to establish sm all inclinations away from edge-on (jKalas 


& Jewitt 1995). A small inclination away from edge-on 
(85-89"") IS dithcult to ascertain at large separations be¬ 
cause the sharpness of the disk midplane in projection 
(i.e., the shape of a cut perpendicular to the midplane) 
is a combination of the disk scale height and the small 
inclination to the line of sight. Closer to the star, how¬ 
ever, the small inclination combined with an asymmetric 
scattering phase function tends to shift the isophotes so 
that the disk does not exactly intersect the star. For 
example, if the disk midplane appears to pass “above” 






























































































































P Pic’s inclined disk in polarized light 
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Table 1 

GPI observations of (3 Pic 


Date 

Observing 

Mode 

Exposure 
Time (s) 

Parallactic 

Rotation (°) Seeing (") 

d Pic b 

Separation (mas) PA (°) 

2013-11-16 

Kl-Spec. 

1789 

26 

1.09 

430.3 ± 3.2 

212.31 ±0.44 

2013-11-16 

K2-Spec. 

1253 

18 

0.93 

426.0 ±3.0 

212.84 ±0.42 

2013-11-18^ 

H-Spec. 

2446 

32 

0.68 

428.1 ±2.7 

212.22 ±0.39 

2013-12-10 

H-Spec. 

1312 

38 

0.77 

418.8 ±3.6 

212.64 ±0.53 

2013-12-10^ 

J-Spec. 

1597 

18 

0.70 

419.1 ±6.2 

212.16 ±0.81 

2013-12-11 

H-Spec. 

556 

64 

0.46 

419.2 ±5.1 

212.26 ±0.72 

2013-12-12 

H-Pol. 

2851 

91 

- 

426.6 ± 7.0 

211.80 ±0.68 

2014-03-23 

Kl-Spec. 

1133 

47 

0.47 

412.5 ±2.7 

212.08 ±0.41 

2014-11-08 

H-Spec. 

2147 

25 

0.77 

362.9 ±4.1 

212.17 ±0.65 

2015-01-24 

H-Spec. 

716 

5 

0.85 

347.7 ±4.7 

212.17 ±0.65 


^ These observations were published by [Macintosh et al.| ( |^14| ), but have been re-reduced here to maintain 
homogenaeity across the datasets _ _ 

^ These observations were published by |Bonnefoy et al.| ( |2014| ), but have been re-reduced here to maintain 
homogenaeity across the datasets 


the star, then that is taken as evidence that the disk 
conies out of the sky plane above the star, and enhanced 
forward scattering leads to the apparent misalignment 
between the midplane and the star. 

For P Pic, |Milli et al.| (|2014|) discovered that the disk 
midplane traces a line that lies above the star. They 
inferred an ^86-89° disk inclination from modeling the 
data, using a Henyey-Greenstein phase function with g 
= 0.36. One significant issue with inferring the line-of- 
sight inclination from their V dataset is that the 3.8 /im 
morphology of the disk within 10 AU is a combination of 
scattered light and thermal emission. Therefore the very 
warm dust n ear the star contrib utes to the detected flux 
within 0.5". Milli et al. (2014) concluded that shorter 


wavelength observations, that are less-contaminated by 
thermal emission, are necessary to disentangle the geom¬ 
etry of the system within 0.5" radius. 

The technique used to image 3 Pic b relies on an- 


gular differential imaging (A 
achieve sub-arcsecond inner 

DI; |Marois et al.l 

20061 ) to 

working angles (f 

Lagrange 

et al. 2012 Milli et al. 

2014 

1 Nielsen et al.| 20141). Pbr 

ground-based observations, tl 

iis technique typically pro- 


vides more effective point spread function (PSF) sub¬ 
traction than using PSF reference stars images, which 
are subject to the time variability of the AO-corrected 
PSFs. However, when applied to extended objects— 
such as circumstellar disks—ADI often c auses significant 
self-subtraction 


(e-g- 


2012 ), impacting the 


Milli et al 

accuracy of derived model disk parameters. These ef- 
fects can be mitig ated with forward modeling (e.g. \Es 


posito et al. 2014), but self-subtraction can be largely 
avoided through polarimetric differential imaging (PDI; 
Kuhn et ar]|2001 ). PDI takes advantage of the fact that 


scattered light is inherently polarized while stellar ra¬ 
diation is not, to subtract the unpolarized stellar PSF, 
revealing the polarized disk underneath. 

Here we present polarimetric observations of P Pic’s 
debris disk at 1.6 /im (i7-band), taken with the Gemini 
Planet Imager (GPI). The data simultaneously reveal the 
debris disk in polarized light and P Pic h in unpolarized 
light. These observations provide a unique perspective 
on the vertical extent of the disk at small angular sepa¬ 
rations, where ADI self-subtraction is typically the most 


severe. In addition, we present new astrometric mea¬ 
surements of the companion P Pic b taken with GPI’s 
spectroscopy mode, which we use to provide an updated 
orbital fit. 

In ^we provide a description of the observations and 
data reduction steps for both polarimetry and spectral 
mode data. We describe our analysis of the disk image 
in ^ and our orbit fitting in ® In ^ we discuss our 
interpretation of the two fits, l^th individually and in 
the context of the disk-planet interaction. We present 
our conclusions in ^ 

2. OBSERVATIONS AND DATA REDUCTION 

GPI is a recently commissioned NIR instrument on the 
Gemini South telescope, designed specifically for t he di¬ 
rect imaging of exo planets and circumstellar disks (Mac¬ 


intosh et al. 


The optical path combi nes high- 


.2014). _ 

order adaptive optics (AO; [Poynee r et a'r]|2QI4|), with an 
apodized pupil Lyot coronagraph ( |Sonmmer et al.||2QII|) 
field spectrograph (H'S; Larkin 


et al. 

[2014) 

centra 

1 star 


The coronagraph system masks out the 
while simultaneously suppressing diffrac¬ 
tion caused by the telescope and its support structure. 
Within the IFS, GPI’s focal plane is sampled by a lenslet 
arr ay at a spatial scale of 14.13 mas/lenslet (see Sec¬ 
tion 4.1) over a ^ 2.8" x 2.8" square field of view. The 
light from each lenslet is passed through either a spec¬ 
tral prism, to allow for low resolution {R ^45) integral 
field spectroscopy, or a Wollaston prism, for broadband 
integral field polarimetry. During observations, Gemini’s 
Gassegrain rotator is turned off to allow the sky to ro¬ 
tate while the orientation of the PSF remains static with 
respect to the instrument. 

The complexity of the instrument results in an intri¬ 
cate path from raw data to a fully processed datacube, 
requiring many calibrations and transformations to ob¬ 
tain a final data product. As a result, the GPI Data 
Reduction Pipeline (DRP) has been designed as a ded¬ 
icated software application for processing GPI d ata. A 
full descript i on of the GP I DRP can be found in Perrin 
et al. (2014), Maire et al. (20I0|) and references therein, 
walkthrough of the data reduc tion process fo r GPI 


polarimetry data can be found in Perrin et al. (2015) 
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Figure 1. Left: The (3 Pic debris disk in polarized intensity, Qr, from observations with GPPs polarimetry mode. The disk image has been 
rotated so that the midplane of the outer disk (P.A.= 29.1°; black dashed horizontal line) is horizontal. The star’s location and /3 Pic 6’s 
location are marked by the magenta x and black circle, respectively. Center: The Qr image with polarization vectors overplotted. Though 
the centrosymmetric nature of the polarization is captured in the transformation to the radial stokes images, the vectors serve to emphasize 
this behavior. Right: The radial polarized intensity, Ur, from the same datacube as the image on the left, shown at the same color scale. 
For optically thin circumstellar material, the flux is expected to be solely in the Qr image, thus this images provides an estimate of the 
noise in the disk image. 


Below, we provide a brief summary of the observations 
and relevant data reduction steps. All the data described 
herein was reduced using the GPI Pipeline version 1.2 or 
lateiQ 

2 .1. Polarimetry mode observations 

Polarimetric observations of p Pic were carried out on 
2013 December 12 UT, while performing a series of AO 
performance and optimization tests (Table [^. [3 Pic was 
observed for a total of forty-nine 60 s frames, during 
which the field rotated in parallactic angle by 91°. Be¬ 
tween each image the half waveplate (HWP) modulator 
was rotated by 22.5°. For 25 frames, GPP s two Ster¬ 
ling cycle cryocoolers (Chilcote et al. 2012|) were set to 
minimal power to reduce vibration in the telescope and 
instrument to improve the AO performance. The exter¬ 
nal Gemini seeing monitors were not operational during 
these observations and as a result the seeing throughout 
the sequence remains unknown. However, ^ Pie h can 
easily be seen in the majority of raw detector images, 
even before data reduction. 

Each raw data frame was dark-subtracted, corrected 
for bad-pixels and then ‘destriped’ to remove any re¬ 
maining correlated noise i n the raw image caused the by 
the cryocooler vibration (Ingraham et al. 2014). Since 
the time of these observations, the level of vibration has 
been significantly mitigated through the use of a new 
controller, which driv es the two coolers 180° out of phase 
(jHartung et al.||2014|). The vibration caused by the two 
coolers now interferes destructively and the overall effect 
is significantly damped. With a reduced level of vibra¬ 
tion, the destriping algorithm is only needed for very 
short exposure times, and incorrect use may result in 
the injeetion of unwanted noise. We direct the reader to 
the GPI IFS Data HandboolEl for further details on the 
destriping algorithm and its appropriate use. 

In GPPs polarimetry mode (pol. mode), a Wollaston 
prism splits the light from each lenslet into two spots 


of orthogonal polarization states on the detector. Flex¬ 
ure effects within the instrument cause these lenslet PSF 
spots to move from their predetermined locations on the 
detector, typically by a fraction of a pixel. For each 
frame, the PSF offset was determined using a cross¬ 
correlation between the raw frame and a set of lenslet 
PSF models measured using a Gemini Facility Calibra¬ 
tion Unit (GC AL) calibration fram e. The overall method 
is decribed in Draper et al. (2014) using high-resolution 


microlenslet PSfl's. Here we use a Gaussian PSF model 
which is less computationally intensive, but provides sim¬ 
ilar results. The raw frames were then reduced to po¬ 
larization datacubes (where the third dimension carries 
the two orthogonal polarizations) using a weighted PSF 
extraction centered on the flexure -corrected location of 
each of the lenslets’ two spots (see Perrin et al. |2Q15 ). 

Each cube was divided by a reduced GCAL flat field 
image, smoothed using a low pass filter. The flat field 
corrects simultaneously for throughput across the field 
and a spatially varying polarization signal. In theory, 
this polarization signal should be removed during the 
double differencing procedure later in the pipeline; how¬ 
ever, we have found empirically that this polarization 
signal is best divided out of each cube individuall}0 

To determine the position of the occu lter-obscured 
star, a Radon-transform-based algorithm (jPueyo et al. 
'20151) was used to measure the p osition of the elongated 


satellite spots (Wang et al. 2014). Knowledge of the ob¬ 
scured star’s location is critical when combining multi¬ 
ple datacubes that must be both registered and rotated. 
Each datacube was then corrected for distortion across 
the field of view (jKonopacky et al.|2014|) . The datacubes 
were then corrected for any non-common path biases be¬ 
tween the two polarization sp ots using the doubl e differ¬ 
encing correction described by 
being smoothed with a 2-pixe! 


Perrin et al. (2015), before 
fl'WHlVl Gaussian profile. 
Instrumental polarization, due to optics upstream of 
the waveplate, converts unpolarized light from the stel- 


^ http://planetimager.org/datapipeline 
^ http://docs.planetimager.org/pipeline/ifs/index.html 


^ This feature will be included in the release of GPI Pipeline 
version 1.4. 
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lar PSF into measurable polarization that, if left uncal¬ 
ibrated, can mimic an astrophysical signal. This signal 
was removed from each difference cube individually, first 
by measuring the average fractional polarization (i.e., the 
difference of the two orthogonal polarization slices di¬ 
vided by their sum) inside of the occulting mask, where 
the flux is due solely to star light diffracting around the 
mask. We assume that this fractional polarization signal 
is due to polarization of unpolarized stellar flux by the 
instrument and telescope. For each lenslet, the fractional 
instrumental polarization was then multiplied by the to¬ 
tal intensity at that location, and then subtracted off in 
a similar manner to the double differencing correctiorj^ 
Using this method we find the instrumental polarization 
to be ^ 0.5%, a sini ilar level to that reported by |Wik-| 


torowicz et al. (2014) using the same dataset. 

'I’he ditterence cubes were then shifted to place the ob¬ 
scured star at the center of the frame and then rotated 
to place north along the ^-axis and east along the x-axis. 
All of the polarization datacubes were then combined 
using singular value decomposition matrix inversion to 
ob tain a three dimen sional Stokes cube, as described 
in Perrin et al.| (|2015 ). Non-ideal retardance in GPFs 
HWP makes GPi weakly sensitive to circular polariza¬ 
tion, Stokes V. Measurements of the circular polariza¬ 
tion of an astrophysical source would require knowledge 
of the HWP’s retardance beyond the current level of cali¬ 
bration. Therefore, in almost all cases the Stokes V cube 
slice should be completely disregarded. 

The Stokes datacube was then transformed to ‘ra- 
dial’ Stokes parame ters: (/, Q, U, U) ^ (/, (5r,Ur,U) 


( Schmid et al.||2QQ6[ lp' Under this convention, each pixel 


in the CJ^ image contains all the linear polarized flux 
that is aligned perpendicular or parallel to the vector 
connecting that pixel to the central star. A positive Qr 
value indicates a perpendicular alignment and a negative 
Qr indicates a parallel alignment. Note that this sign 
convention is opposite that used in|Schmid et al.|(|2006|), 
where positive values of Qr correspond to a parallel align- 
ment. The Ur image holds the flux that is aligned ±45° 
to the same vector. For an optically thin circumstellar 
disk, the polarization is expected to be perpendicular and 
all the flux is expected to be positive in the Qr image. 
The Ur image should contain no polarized flux from the 
disk and can be treated as a noise map for the Qr image. 
The final reduced disk image can be seen in Figure 

2 .2. Spectral mode observations 

Observations of (3 Pic in spectroscopic mode (spec, 
mode) were carried out during four separate GPI com¬ 
missioning runs, as well as during an ongoing astrometric 
monitoring program scheduled during regular general ob¬ 
serving time. In total, we present nine individual sets of 
observations over seven epochs (Table f^. Two of the 
observation sets have been previously p ublished: the H- 
band dataset from 2013 November 18 (Macintosh et al. 
2Q14|) and the J-band dataset from 2013 December 10 
( [Bonnefoy et al.|[~2014 ). Here we have re-reduced all the 
data in a consistent manner in an effort to reduce system¬ 
atic biases and maximize the homogeneity of the dataset. 
As with the polarization mode observations, those obser¬ 
vations that were taken during the instrument’s commis¬ 
sioning were carried out during AO optimization tests 
and therefore have a range of exposure times and filter 


combinations. 

All datasets were reduced with standard recipes pro¬ 
vided by the GPI DRP. Raw data frames were dark sub¬ 
tracted and destriped for microphonics in the same man¬ 
ner as the polarimetry observations. A short-exposure 
arc lamp image was taken contemporaneously with each 
science observation to measure the offsets of the lenslet 
spectra due to flexure within the IFS. The mean shift was 
calculated for a subset of lenslets across the field of view 
relative to a high SNR arc lamp image taken at zenith 
via a Leve nberg-Marquardt least-squares minimization 
algorithm (jWolff et al.||2014|). 

The raw detector iniage was then transformed into a 
spectral datacube, using a box extraction method. For 
observations obtained with the K1 and K2 filters, ther¬ 
mal sky observations were taken immediately before or 
after the observation sequence. Sky background cubes 
were created in the same manner described above and 
subtracted from science dat acubes. Finally, all cu bes 
were corrected for distortion (jKonopacky et al.||2014). 

Each dat a-set was PSF subt racted using the methods 
outlined in Pueyo et al. (2015). To minimize systematic 
biases, the ensemble of data-sets was treated as uniformly 
as possible. The main steps of this data reduction process 
include: high-pass filtering, to remove the remaining PSF 
halo; wavelength-to-wavelength and cube-to-cube image 
registration, to correct for atmospheric differential re¬ 
fraction and sub-pixel stellar motion across the observ¬ 
ing sequence; subtracting the speckles using the KLIP 
princ ipal component analysis algorithm (Soummer et al 


2012 ) on each wavelength slice in each cube; rotation to 
align the north angle of each image; and co-adding the 
resulting cubes in time. 

For the epochs in which P Pic b was observed on con¬ 
secutive nights, relative alignment was test ed using both 
the cr oss correlation method described in |Pueyo et al.| 
(2015) and the absolute stellar locations based on the 
satellite spot centroids derived using the GPI DRP. For 
these epochs we found better consistency in the location 
of Pic b using the DRP centroids, which we then chose 
to adopt for all datasets. 

The KLIP algorithm was impl emented using both 
spectral differential imaging (SDI; jMarois et al. 2000) 
and ADI, building for each slice a PS4' library that takes 
advantage of the radial and azimuthal speckle diversity 
(in wavelength and in PA, respectively). Due to the rela¬ 
tive brightness of /3 Pic b with respect to the neighboring 
speckles we limited the exploration of KLIP parameter 
space to two zone geometries and two exclusion criteria 
(1 and 1.5 PSF FWHM) for each dataset. For each slice, 
the 30 PSFs that were the most correlated in the region 
where /3 Pic b is located were used for PSF subtraction, 
except for the J-band data which required 50 PSFs for 
satisfactory subtraction. 

To determine the optimal number of principal com¬ 
ponents to use for each dataset, we examined both the 
evolution of the extracted spectrum and the astromet¬ 
ric stability as a function of wavelength as we increased 
the number of components. This latter test helps us to 
rule out the pathological cases for which either a resid¬ 
ual speckle (i.e. insufficiently aggressive PSF subtrac¬ 
tion) or self-subtraction (i.e. over-aggressive PSF sub¬ 
traction) bias planet centroid estimates. We checked 
for potential biases by comparing astrometric positions 
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Figure 2. Left: The residuals of the Qr image after subtrac¬ 
tion of the best fit disk model, shown at a stretch to emphasize 
their structure. The residual structure in the NW-SE direction is 
likely due imperfect subtraction of the instrumental polarization, 
but may also represent the difference between the true scattering 
properties of the grains and the Henyey-Greenstein function used 
in the model. Right: The Ur image, shown at the same stretch as 
the residuals. The star’s location is marked with a magenta x. 


measured when using only a high-pass filter with those 
measured when applying KLIP. Finally, we checked for 
self-consistency by injecting six synthetic point sources 
at the same separation as /3 Pic 6, but at different posi¬ 
tion angles. Based on these tests we concluded that the 
astrometric measurements do not feature systematics ei¬ 
ther introduced by residual speckles or biases associated 
with KL IP above the uncertainty levels reported in Sec¬ 
tion [FTI 


3. DISK RESULTS 


The debris disk is recovered in polarized light from 
^ 1.7" (32 AU), to an inner working angle of ^ 0.32" 
(6.4 AU); see FigureWhile GPFs iL-band focal plane 
mask extends to a radius of ^ 0.12", uncorrected in¬ 
strumental polarization and other noise sources domi¬ 
nate over the disk emission at separations smaller than 
- 0.32". 


A comparison of the Qr and Ur images indicate that 
the disk is detected at a high signal-to-noise ratio (SNR) 
out to the edge of the GPI field. Overplotting linear po¬ 
larization vectors indicates that the emission is linearly 
polarized perpendicular to the scattering plane, as ex¬ 
pected for optically thin conditions. This property is 
captured in the transformation to radial Stokes param¬ 
eters, but we have included the vectors in Figure for 
additional clarity. 

Morphologically, the disk appears vertically offset from 
the midplane of the outer disk in the NW direction, in¬ 
dicative of a slight inclination relative to the line of sight. 
This is consistent with previou s models of the disk at 
similar angular separations (e.g. Milli et ar]|2Q14 ). 

The Ur image shows low level structure in the form of 
a dipole-like pattern with positive emission in the E-W 
direction. Figure displays the Ur image with a color 
scale that emphasizes this structure. In the radial Stokes 
basis, this is the pattern produced by a constant linear 
polarization across the field, which could be associated 
with residual instrumental polarization that was not suc¬ 
cessfully subtracted during the data reduction process. 
Since the level of these residuals is much lower than the 
disk emission, we defer improvement of our instrumental 
polarization subtraction procedure for future work. 

The disk is not detected in total intensity (Stokes /; 
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Figure 3. Top Left: The Stokes / image from the polarimetry 
observations, without PSF subtraction. The dashed line indicates 
the position angle of the outer disk. The planet can be seen at a 
separation of ~ 0.4" just above the horizontal line, to the SW from 
the central star. Top Right: The Stokes / image after applying 
KLIP/ADI PSF subtraction. The planet is recovered at a very high 
SNR. Bottom Left: The polarized intensity image, Qr, after disk 
model subtraction. The black circle indicates the location of the 
planet in the Stokes / images. Bottom Right: The radial Stokes 
image (same as in Figure [R, Ur- No point source polarization 
signal is detected for /3 Pic o in either Stokes Qr or Ur- All images 
have been rotated so that the outer disk’s PA is horizontal (dashed 
black line). In all four images the star’s location is marked with a 
magenta x. 

Figure]^, where images are dominated by the residual 
uncorrected PSF of the star itself. Due to both the ex¬ 
tended nature of the disk at these angular scales and 
frame-to-frame variation of the PSF (compounded by 
the AO tests carried out during the observing sequence), 
ADI PSF subtraction has proven unsuccessful. Without 
an unbiased total intensity image of the disk, characteri¬ 
zation of the polarization fraction remains out of reach at 
present. As a result, we opt to model only the polarized 
intensity. 


3.1. Disk modeling 

The principal objective of our disk modeling is to re¬ 
trieve basic geometric properties of the disk. The model¬ 
ing approach adopted here is to combine a simple recipe 
for the 3D dust density distribution with a parametric 
model of the polarized scattering phase function and then 
fit to the data using the affine-invariant sampler in a 
parallel-tempering scheme from the emcee Markov chain 
Mont e Garlo (MGMG) package (Foreman-Mackey et al. 
2013). Parallel tempering uses walkers at difterent Tern 


peratures’ to broadly sample the posterior distributions 
and is therefore a useful strategy when the likelihood 
surface is complex. 

For a disk seen in edge-on projection, the radial dust 
density distribution becomes degenerate with the dust 
scattering properties. This degeneracy is typically bro¬ 
ken with the use of physical grain models, which describe 
scattering properties (including polarization) as a func- 
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Table 2 

Disk Model Parameters 


Parameter 

Symbol 

Range 

Prior Distribution 

Inner Radius 

Ri 

1 - lOOAU 

Uniform in logi^i 

Outer Radius 

R 2 

Ri - 500AU 

Uniform in logi^2 

Density Power Law Index 

h 

0.5-4 

Uniform in (3 

Scale height aspect ratio 

ho 

0.01-2 

Uniform in log ho 

HG asymmetry parameter 

9 

0-1 

Uniform in g 

Line of Sight Inclination 

0 

80 - 90° 

Uniform in coscp 

Position Angle 

OpA 

25 - 35° 

Uniform in Op a 

Flux Normalization 

No 

1-1000 

Uniform in logA^o 

Flux Offset 

lo 

-5 - 10 

Uniform in /q 


tion of wavelength. In practice, observations are fit to 
grain models either using simul taneous polarization and 
total intensity info rmation (e.g. Graham et ah 2007), or 
multicolor images (|Golimowski et al.||2006|). With only 
single wavelength polarized intensity images available, we 
inste ad use the Henyey-Greenstein (HG) scattering func¬ 
tion (Henyey & Greenstein 1941) to describe the scat¬ 
tering eHiHenc}71[i"iriimHi(m^ scattering angle. The 
shape of the HG scattering function is a function of only 
one parameter, the expectation value of the cosine of the 
scattering angle, g = (cos 0 ), and thus provides a useful 
tool to approximate grain scattering when using physical 
models is impractical. The applicability of the HG scat¬ 
tering function to the modelin g of our polarized intensity 
images is discussed in Section [5T| 

Our dust density model, expressed in stellocentric co¬ 
ordinates, 77(r, 2 ;), follows a power law between an inner 
radius, i?i, and an outer radius, i? 2 , and has a Gaussian 
vertical profile with RMS height h^r and constant aspect 
ratio, ho'. 


? 7 (r, z) oc 





where r is the radial distance from the star, z is the 
height above the disk midplane and P is the power law 
index of the dust density. Inside Ri and outside R 2 the 
dust density is zero. The dust density distribution is 
combined with the Henyey-Greenstein function, 
to generate a scattered light image of the disk as seen 
in 2D projection from the observer’s frame, where the 
intensity for a given pixel {x'^y') is calculated as the 
integral along the line-of-sight direction z'‘. 

= I 0 + [ ^ri^{r,z)H{0,g)dz'. 

Jz' = -R2 ^ 

Here, 7^0(r, z), represents the dust density distribution, 
but tilted with a disk inclination, 0, relative to the ob¬ 
server’s line of sight. The scattering angle 0 = 9{x'^y'^ z') 
is a function of position. 

The 1/r^ term accounts for the diminishing stellar flux 
as a function of distance from the star. The disk’s po¬ 
sition angle, Opa^ is implemented as a coordinate trans¬ 
formation between the stellocentric coordinates and the 
projected observer’s coordinates. The constant, Jq, and 
the flux normalization, Nq, have been included to ac¬ 
count for any possible biases and the conversion between 
model flux and detector counts, respectively. In sum¬ 
mary, our model has a total of nine free parameters: 


Ri, R 2 , /^, ho,g, 0, OpA, No, Jq. 

Within the current model there exists a degeneracy be¬ 
tween forward scattering (^ > 0) with an inclination of 
< 90° and backwards scattering {g < 0) with an incli¬ 
nation of (j) > 90°. In an effort to conserve computation 
time we chose to assume forward scattering and place 


parameter, q > 0 , 


a prior constraint on the scattering 
which is consistent with the model of Milli et al. (|2014|). 
A summary of the model parameters and their prior dis- 
tributions can be found in Table [21 

We fit the model to the GPI disk image using the 
parallel-tempering sampler from the emcee package. The 
diffraction limit in the i7-band for Gemini south is 
0.043", equal to about three GPI pixels. We therefore 
apply a 3x3 pixel binning to both the Qr and Ur im¬ 
ages before fitting. This improves the noise statistics 
and speeds up the execution time of the MCMG fit, with¬ 
out significantly sacrificing spatial information. At each 
angular separation in the Qr image, the errors were esti¬ 
mated as the standard deviation of a 3 pixel-wide annulus 
centered at that separation in the Ur image. The error 
estimates therefore contained photon noise, read noise 
and the unsubtracted instrumental polarization. 

The MCMG sampler was run for 2500 steps with 2 
temperatures, 128 walkers and burn-in of 500 steps. Ad¬ 
ditional temperature chains were not employed because 
of the additional computational cost incurred and the 
lack of evidence that the Markov chain sampler was only 
selecting local islands of high likelihood. One strength 
of using ensemble sampling over other types of sampling 
for MCMG fitting is that large speed-ups are possible via 
parallel-processing. On a 32-core (2.3 GHz) computer 
the entire MCMG run took nearly five days to complete. 

After the run, the maximum auto-correlation across 
all parameters was found to be 85 steps, indic ating that 


the chains should have reached equilibrium (Foreman- 
Mackey et al.|[2M3 recommend ^ 0(10) autocorrelation 


times for convergence). In addition, the chains were ex¬ 
amined by eye and appeared to have reached steady-state 
by the end of the burn-in phase. The posterior distribu¬ 
tions (Figure were estimated from the zero temper¬ 
ature walkers, using only one of out every 85 steps to 
ensure statistical independence of the surviving samples. 
T he expected covarianc e between the inclination 0 and 
g (jKalas & Jewitt|1995|) is reproduced. Degeneracies are 
also found between 0 and g, ho and g, ho and [3, ho and 
(j), /3 and g, and [3 and 0, but in all cases, the parameters 
appear to be well constrained. 

The 16%, 50% and 84% percentiles for each parameter 
are displayed in a table in the upper right corner of Fig- 
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Parameter 

Percentiles 

16 

50 

84 

Ri[AU] 

23.04 

23.61 

24.47 

R2{AU] 

126.18 

138.84 

158.12 

/? 

0.71 

0.85 

1.01 

ho 

0.131 

0.137 

0.142 

9 

0.729 

0.736 

0.744 


85.08 

85.27 

85.53 

dpA[°] 

30.22 

30.35 

30.49 


T r 



Figure 4. The posterior distributions of the model parameters from MCMC disk model fitting to the Qr disk image. The diagonal 
histograms show the posterior distributions of each parameter marginalized across all the other parameters. In each plot, the dashed lines 
indicate the 16%, 50% and 84% percentiles. The off-diagonal plots display the joint probability distributions with contour levels at the 
same percentiles. The normalization term, Nq, and the constant offset term, /q, have been excluded from this plot because they convey no 
relevant astrophysical information. Inset table: The 16%, 50% and 84% percentiles from the marginalized distributions. 


ure[^ Marginalized across all parameters we find a disk 
inclined relative to the line of sight by 0 = 85.27 °Iq! 19 ? 
with an inner radius of Ri = 23.6lo!6 an outer 

radius of R 2 = AU and an aspect ratio of 

ho = 0.137lo‘.oo6* position angle of the disk is 

OpA = 30.35°to 28 5 where the errors include GPI’s sys¬ 
tematic error in position angle 0.2°). Note that the 
systematic uncertainty in the position angle is not re¬ 
flected in Figure The dust is well fit by forward scat¬ 


tering grains, with a scattering asymmetry parameter of 
g = 0.736lo.oo7- These results are further discussed in 
Section 5.1. 

Figure [^displays the best fit model and the residuals 
of the Qr image minus the model. The best fit model 
was generated using the median value of each parameter 
in the marginalized posterior distribution. We find that 
the highest likelihood disk model successfully reproduces 
the GPI data. When examined at a different color scale, 
the residuals image displays similar low-level structure 
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Figure 5. Left: The disk model generated with the median values from the marginalized posterior distributions (as found in Figure 
The inner edge of the disk is at a projected separation of 1.2", but contributes negligible light to the observed surface brightness. Cen\ 
The residuals of Qr image minus the disk model. The level of the residuals is very similar to the Ur image. Right: The Ur image, 
reproduced here as a point of comparison to the residual image. In all images the star’s location is marked with a magenta X .The images 
have been rotated so that the outer disk is horizontal and all are displayed at the same colour scale as Figure^ 

as that of the Ur image (Figure]^. The structure in the 
NW-SE direction is likely the Qr counterpart of the resid¬ 


ual instrumental polarization that’s seen in the Ur image. 
A possible alternative explanation is that the structure 
could be due to a mismatch between the true scattering 
properties of the dust and the Henyey-Greenstein scat¬ 
tering function at small angular separations. A second 
structure can be seen along the disk midplane to the NE 
of the star. This asymmetric brightness feature is possi¬ 
bly due to a local overdensity of dust, that would increase 
the scattering at that location. Indeed, the /3 Pi c disk 


is known to have multiple brightness asymmetries (Apai 
et al.j[2M5 provide a good summary). However, the tW 


ture IS detected at similar brightness levels as the resid¬ 
ual instrumental polarization and may yet be an unchar¬ 
acterized artifact of the polarimetry reduction. Deeper 
observations of the disk will be required to distinguish 
between a true brightness asymmetry and instrumental 
effects. 

4. PLANET RESULTS 
4.1. Astrometry in Spectroscopy Mode 

We describe here in broad terms our astrometric mea¬ 
surements and estimation of uncertainties, without delv¬ 
ing into the details of each individual dataset. Eor each 
epoch, the entirety of the dataset is combined to estimate 
the planet’s position relative to P Pic. The errors on this 
relative position are a combination of the error on the 
star’s position, the planet’s position, GPI’s pixel scale 
and the accuracy to which we know GPI’s orientation 
relative to true North. 

Eor each dataset, the stellar position was calculated 
using two methods. The wavelength slices of each dat- 
acube were first registered usin g th e relative alignment 
procedure described in Section |2.2| and then collapsed 
into a broadband image. A Radon transform was then 
performed on the radially elo ngated satellite sp ots to 
find the stellar position (as in Pueyo et al. 2015). The 
stellar position was also estimated using the geometric 
mean of the satellite spot locations provided by the GPI 
DRP. Most i^-band datasets show agreement between 


two methods at the 0.05 pixel level, with the exception of 
the 2013 Dec. 11 commissioning sequence, during which 
extensive AO performance tests where being carried out. 
Eor iU-band data-sets the difference between the two 
methods is no more than 0.05 pixels and for the J-band 
data-set it is 0.2 pixels. We found greater consistency 
in the relative location of p Pic b between observations 
obtained on consecutive nights when using the Radon 
method, and therefore chose to adopt the centroids mea¬ 
sured with the Radon method for all measurements. Eor 
each dataset, we considered the difference between the 
two methods as our estimate for the uncertainty on stel¬ 
lar position. 

The location of P Pic h (in detector coordinates) was 
estimated at each wavelength channel, at each epoch and 
in each filter using the modified matched filter described 


Pueyo et al. (2015). The uncertainty in P Pic 6’s 


location was estimated as the scatter in the position of 
the planet as a function of wavelength and number of 
principal components. We found the uncertainly to range 
from 0.05 pixels, for the datasets with significant field 
rotation and where the planet was at larger separations, 
up to 0.15 pixels, for the later epochs where the planet 
is significantly closer to the stellar host and SDI is less 
effective. 

We esti mated GPI’s pixel scale u sing the methods de¬ 
scribed in Konopacky et al. (2014) by combining all the 
data presented therein with four new observations of 0^ 
Ori B, taken between September 2014 and January 2015. 
We find an updated pixel scale value of 14.13 ± 0.01 
mas/lenslet. Konopacky et al. (2014) measured a PA 
offset of —1.00 ± 0.03"" during GPl cornmissioning. Sub¬ 
sequently, version 1.2 of the GPI DRP was updated to 
incorporate that 1° offset and correct for it automati¬ 
cally. Using the new measurements of 0^ Ori B, we find 
a residual PA offset of —0.11 ± 0.25°. 

Based on the measured location of p Pic b and its par¬ 
ent star in detector coordinates we calculated the relative 
separation and position angle at each epoch. The sep¬ 
aration was converted to milliarcseconds using the new 
platescale estimate and the PA was adjusted by —0.11°. 




















10 


Millar-Blanchaer et al. 


The separation and position angle from each measure¬ 
ment can be found in Table [H Uncertainties on these 
quantities were combined with the errors on the star po¬ 
sition and planet position to yield the errors presented 
in the table. 

4.2. Astrometry in Polarimetry Mode 

P Pic b is detected in the Stokes / image as a point 
source superimposed on the extended PSF halo (Fig¬ 
ure |^. After applying PSF su btraction using a python 
impfenentation of KLIP/ADI (|Wang et al.||2015| to the 
image, the planet is recovered at extremely high SNR. 
The planet’s position in the Stokes / image was estimated 
using the StarFinder IDL package (Diolaiti et al.|2000), 
which requires the user to input a PSF' model for preci- 
sion astrometry. In GPFs polarimetry mode the entire 
bandpass is seen by each frame and therefore the satel¬ 
lite spots are elongated and cannot be used as a PSF 
reference, as they are in spectroscopy mode. Instead, we 
used a GPI PSF generated w ith AO simulation software 
( Poyneer fc Macintosh||2QQ6 ). 

I'o estimate astrometric errors we used StarFinder to 
measure (3 Pic 6’s location in the total intensity image 
from each of the 49 polarization data cubes. The RMS 
difference between the planet location in the individual 
cubes and the Stokes / image was taken to be the error 
in the planet location. The error on the location of the 
star is estimated from the RMS scatter of the measured 
star’s position across the set of cubes. This error tracks 
the motion of the star behind the coronagraph between 
frames, which we expect to be larger than the errors on 
the star’s position determined by the Radon transform, 
and therefore likely overestimates the errors. 

The position of p Pic b in the polarimetry mode ob¬ 
servations can be found in Table As with the spec¬ 
troscopy mode data, the errors represent a combination 
of the errors on the star’s and planet’s positions, GPFs 
pixel scale and GPFs PA offset on the sky. 

4.3. Orbit fitting 

Using the ten newly obtained astrometric points (nine 
from spec, mode and one from pol. mode), combined 
with the datasets presen ted by Ghanvin et al.| ( |2QI2 ) 


and Nielsen et al. 
orbital elements o 


(2014), we fit for the six Keplerian 
P Pic b plus the total mass of the 


system using the parallel-teni pered sampler from emcee 
(Foreman-Mackey et al.||20I4). While astrometric data- 
points have been published in other papers, in an effort to 
minimize systematics between datasets, we limited our¬ 
selves to only these two large datasets where considerable 
effort has been made to calibrate ast rometric erro r s. Th e 
fitting code was previo usly used in iKalas et ah (2013), 


Macintosh et al. (2014), and Pueyo et al. 


ht the radial velocit y measurement of 


VVealso 
;ne planet from 


Snellen et al. (2014), which allows us to constrain the 
line-ol-sight orbital direction and break the degeneracy 
between the locations of the ascending and descending 
node. 

The model fits seven parameters: the semi-major axis, 
a; the epoch of periastron, r; the argument of periastron, 
00 ; the position angle of the ascending node, U; the in¬ 
clination, i; the eccentricity, e; and the total mass of the 
system, Our orbital frame of reference followed the 


binary star sign convention used in Green (1985). Under 
this convention the ascending node is defined as the loca¬ 
tion in the orbit where the planet crosses the plane of the 
sky (centered on the star), moving southward in projec¬ 
tion. N ote that this is 1 8 0° diff erent from the convention 
used in Ghauvin et al.| (|20I2 ). The projected position 
angle of the ascending node on the sky is defined from 
North, increasing to the East. The argument of the pe¬ 
riastron is defined as the angle between the ascending 
node and the location of the periastron in the orbit, with 
00 increasing from the ascending node. The epoch of peri¬ 
astron, r, is defined in units of orbital period, from 1995 
October 10 (Julian date 2450000.5). A summary of the 
orbital parameters and their prior distributions can be 
found in Tabled 

The MCMC sampler was run for 10,000 steps with 
10 temperatures and 256 walkers after a “burn-in” of 
2000 steps. After the run, the maximum auto-correlation 
across all parameters was found to be 25 steps, indicat¬ 
ing that the chains should have reached equilibrium. The 
posterior distributions (Figure were estimated using 
the zero temperature walkers, using only one of out ev¬ 
ery 25 steps. In Figure the epoch of periastron was 
wrapped around to be only positive between 0 and I and 
the arguments of the periastron was wrapped around to 
range from 0 to 360°. A random selection of 500 accepted 
orbits are plotted on top of the astrometric and radial ve¬ 
locity datapoints in Figures and respectively. While 
the orbital fits are generally consistent with most of the 
astrometric datapoints, the majority of the orbital solu¬ 
tions fall more than I-cr from the measured radial veloc¬ 
ity. 

We find that the planet has a semi-major axis of 
9.25^0 45 AU, an inclination of 89.01° ± 0.30 and an as¬ 
cending node at a position angle of 31.75° ± 0.15. We 
take the median of the marginalized posterior distribu¬ 
tion to be the best estimator of each parameter’s value, 
and the 68% confidence values as the errors. Following 
this convention, the eccentricity of the orbit is found to 
be e = 0.07tg;j5. However, the eccentricity is a positive 
definite quantity and typical estimators (e.g., the mean 
and median) will overestimate the true eccentricity when 
it is small (e < O.I). When considering eccent ricities of 
radial velocity planets, Zakamska et al. (2011) consider 
several different estimators and find that for small ec¬ 
centricities the mode of the distribution is the least bi¬ 
ased. The mode of our distribution falls in the smaller 
eccentricity bin indicating an eccentricity very close to 
zero. Therefore, it is perhaps more appropriate to quote 
the upper limit on the eccentricity e < 0.26 (95% confi¬ 
dence) . 

For orbits with higher eccentricities (e > O.I), the 
epoch and argument of periastron have strong peaks at 
^0.5 periods and ^ 170°, respectively. At lower eccen¬ 
tricities these two parameters remain degenerate, with a 
large range of acceptable values. Overall, the marginal¬ 
ized distributions reveal that these parameters are still 
relatively unconstrained. The ensemble of accepted or¬ 
bits at the end of the run have a reduced of f 

Marginalized across all parameters, the total mass of 
the system is I.6I±0.05Mq. At ^ II Mj^p, Pic b con¬ 
tributes less than 1% to the total mass, giving Pic itself 
a mass of = 1.60 ± 0.5Mq. This falls slightly below 
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Table 3 

Orbit Model Parameters 


Parameter 

Symbol 

Range 

Prior Distribution 

Semi-major axis 

a 

4 - 40AU 

Uniform in log a 

Epoch of Periastron 

r 

-1.0- 1.0 

Uniform in r 

Argument of Periastron 

(jj 

—27r — 27r rad 

Uniform in uj 

Position Angle of the Ascending Node 

Q 

25 - 85° 

Uniform in H 

Inclination 

i 

81 - 99° 

Uniform in cosi 

Eccentricity 

e 

0.00001 - 0.99 

Uniform in e 

Total Mass 

M'j' 

0-3Mo 

Uniform in Mq 


the range estimated by Crifo et al. | ([TM7|) (1.7-1 


and just within the range of 


'Blondel & Djie 
(1.65 — 1.87Mq), who both use evolutionary mode 
the HR diagram to date (3 Pic. 


(2006) 


s and 

This estimate provides 
a slightly sm aller value than that presented in |Nielsen| 
let al.l (|2014|) , though still consistent within their errors 

By combining the semi-major axis and stellar mass val¬ 
ues of each walker at each accepted iteration, we are able 
to create a posterior distribution for the orbital period, 
from which we derive that P = 22.4l^J yr. The large 
upper limit is due to the extended tail in the semi-major 
axis distribution. 


4.4. Planet polarization 
Giant exoplanets may have polarized emission i n the 


ley & Seng 

upta 

Kok et al. 

20T1 


Tor p Pic 6, the recently measur^ 
of ^ 8 hours would induce a polariza¬ 
tion signature due to oblatenes s of less than 0.1% (below 
GPTs current sensitivity limit, Wiktorowicz et al.|2014). 
Therefore, any detected polarization signal would be in- 
dicative of cloudy structure. 

To estimate (3 Pic 6’s polarization, we first created 
a disk-free linear polarized intensity image by combin¬ 
ing the model-subtracted Qr image with the Ur image 
(P = y^Qr 3- The total polarized flux at the loca¬ 
tion of P Pic 6, within an aperture of radius 1.22A/P, was 
then compared to the flux of 38 independent apertures at 
the same angular separation. We find that (3 Pic Ps po¬ 
larized flux is 0 . 5(7 from the mean flux of the independent 
apertures, consistent with zero linear polarization sig¬ 
nal from the planet (see Figure |^. While this measure¬ 
ment does not provide any evidence for cloud structure, 
it does not exclude the possibility either; the magnitude 
of cloud-induced polarization depends on many factors, 
including the atmospheric temperature and pressure pro¬ 
file, the composition, the nature of the inhomogeneities, 
rotation, and viewing angle. The PSF variability dur¬ 
ing the observations makes accurate recovery of the total 
intensity of the planet difficult, and thus we leave the 
characterization of an upper limit on the planet’s polar¬ 
ization fraction for future work. 

5. DISCUSSION 
5.1. The debris disk 

With GPI we probe the projected disk between 0.3" 
and 1.5" at high spatial resolution. The work presented 
here has two advantages over previous attempts to model 
the disk at similar angular separations. First, the polar¬ 


ized intensity images provide a unique view of the disk, 
in particular the vertical extent is free of any biases asso¬ 
ciated with ADI PSF subtraction. Second, the MCMC 
fitting allows us to fully explore the multi-dimensional 
parameter-space and place quantitative confidence inter¬ 
vals on the model parameters. 

MGMG fitting requires evaluation of the likelihood 
function for each set of parameters that is examined. The 
cost of fitting depends on the computational expense of 
evaluating the model and the dimensionality of the model 
parameter space. For that reason we have limited our ex¬ 
ploration to optically-thin scattering, an analytic recipe 
for the phase function, and a simple model of the dust 
distribution. We do not consider mu lti-component disks 


(as modeled for the outer disk, e.g., Ahmic et al. 2009) 
and we assume that the disk aspect ratio is constant. Re- 
gardless of these simplifications, we find that this model 
provides an excellent fit to our polarized image. 

The Henyey-Greenstein scattering function is often 
used to model the total intensity scattering efficiency of 
dust grains, but has not been used extensively for polar¬ 
ized intensity. This is at least partially due to the fact 
that in most circumstances where polarized intensity is 
measured, total intensity is obtained as well, allowing 
for more sophisticated modeling of the dust scattering. 
In addition, the scattering efficiency of polarized inten¬ 
sity of small spherical particles approaches zero at very 
small scattering angles, a feature that is not captured 
by the HG function. While the exact shape of the HG 
function cannot fully reproduce the polarized scattering 
efficiency function for physical models, a quick informal 
survey of possible grain models indicates that our best 
fit = {cosO) « 0.7 can be reproduced by Mie scattering 
particles with a radius of ^ Ijam and an index of refrac¬ 
tion of m = 1.033 —Q.Qlu s i milar to the porous, icy grains 


inferred by|Graham et al. (2007) for AU Mic. However, 
as previously mentioned, a true characterization of the 
physical grain scattering properties will require either an 
unbiased total intensity image, or polarized intensity im¬ 
ages at other wavelengths. We leave the characterization 
of the dust properties of the inner disk f or future work. 

(2014) have a field of 


The observations of Milli et al. 


view (0.4" — 3.8") that overlaps with our disk detec¬ 
tion and therefore provide a good point of comparison. 
They model the L' emission with a single component 
disk model similar to ours, albeit with different radial 
and vertical dust density profiles. Even so, their best fit 
inclination {i = 86°) and position angle {Op a = 30.8°) 
agree fairly well with our own. Their dataset constrains 
the sky-plane inclination less precisely and inclinations 
of 85-88° provide good fits to their data. The consis- 
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a[AU] 


Parameter 

Percentiles 

16 

50 

84 

a[AU] 

8.80 

9.25 

10.71 

r 

0.46 

0.58 

0.69 


75.09 

155.85 

188.92 

141 

31.60 

31.75 

31.90 

41 

88.71 

89.01 

89.31 

e 

0.02 

0.07 

0.18 

Mj 

1.56 

1.61 

1.66 

Period 

20.72 

22.36 

27.69 


Figure 6. The posterior distributio ns of the model para meters from the MCMC orbit model fitting to the astrometry data points, as 
well as the radial velocity value from|Snellen et al.| (|2014|). The diagonal histograms show the posterior distributions of each parameter 
marginalized across all the other pararueters. in eacE plot, the red dashed lines indicate the 16%, 50% and 84% percentiles. The off-diagonal 
plots display the joint probability distributions with contour levels at the same percentiles. Inset table: The 16%, 50%, and 84% percentiles 
from the marginalized distributions. 
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Figure 7. The RA (blue) and DEC (red) offsets of (3 Pic b from 
/3 Pic for a random selection of 100 accepted orbits (dotted lines) 
from the MCMC run. The 29 data points used in the fit are over¬ 
plotted, with the colors indicating their source. Error bars on the 
datapoint are smaller than the mark ers, except for the 2003 mea¬ 
surement fromiChauvin et al.|(|2012 |. This fit includes the_radial 
velocity constraint trom [snelleii et al.| ( |2014| ) shown in Fig. 



Figure 8. Predicted radial velocity of ^ Pic b from from the 
orbital fit constrained by the astrometric data plotted in Fig . fTland 
the single radial velocity measurement from |Snellen et al.| |3oT4|) 
(point with error bar). A random selection ot iOO accepted orbits 
(purple dotted lines) from the MCMC are shown. 

tency between their measurements and ours builds con¬ 
fidence that the measured angles are not highly sensitive 
to the assumed scattering properties and radial dust dis¬ 
tribution. The position angle of the disk seen in our 
images {Opa = 30.35 


-O+0.14 


-0.13. 


and those of 
yned from bot. 

Milli et al. 

1 the outer 

[Apai et al. j 

2015) and the 


mam disk {Opa = 29.1 ± 

warp {OpA ~ 32 — 33°). This ottset, ana now our aisK 
images fit into the context of th e whole system, will be 
further discussed in Section 1531 
The results of our model fitting reveal an inclined disk 
with an inner radius of Ri ~ 23.5 AU (1.2"), populated 
by grains that preferentially forward scatter polarized 
light. The majority of the detected polarized flux is 
therefore inside the projected inner radius and the re¬ 
sult of forward scattering by the constituent dust grains. 


Without direct detections of either the inner or outer ra¬ 
dius, the constraints on both are governed by the overall 
shape and spacing of isophotal contours (see Figure |^. 

The location of the inner edge of the disk seen in our 
model is a unique feature of this work and has not been 
found in previous scattered light imaging at similar an¬ 
gular separations. This could be attributed to both the 
scattering properties of the dust, which make the inner 
edge difficult to see, and the mo deling strategies used in 
those studies. Milli et al. (2014) also use a HG function 
to model their dust. However, their model considers a 
population of parent bodies between 50 AU and 120 AU, 
with the density falling as separate power laws inside 
and outside of these radii and they do not define an in- 
ner radius in the same manner as in our model. |Apai| 


et al. (2015) make surface brightness measurements of 
the disk between 0.5" and 15.0", but find no noticeable 
change in the brightness profile at 1.1". In our model, 
we find that the forward scattering nature of the dust 
grains means that the inner edge itself contributes mini¬ 
mally to the observed surface brightness at its projected 
separation. This serves to emphasize the critical role of 
dust scattering when interpreting the surface brightness 
as a function of radius; a smooth surface density by itself 
does not necessarily exclude features in the radial dust 
profile. 

Note that our model has been defined with a sharp cut¬ 
off inside the inner radius, and caution should be used 
when interpreting the exact value. There may be dust 
inside of the inner radius with a lower surface density. 
For example, the true dust density inside the inner radius 
may have a slowly decreasing i nner power-law, such as 


those considered in Milli et al. (|2014|). 

Imaging and spectroscopic studies in the mid-IR have 
probed similar regions of the debris disk at wavelengths 
where contrast between the stellar flux and the dust 
(thermal) e mission is more favorab le than in the opti¬ 
cal and NIR. Okamoto et al. (2004) found spect roscopic 
evidence for dust belts at b, 10, and 30 AU. |Wahhaj| 


et al. (2003) fit a series of four dust belts to deconvolved 
18 /iui images and found their best fit radii to be 14, 28, 
52, and 82 AU. With the exception of the belts close to 
^15 AU, all of these belts are either well outside or at 
the very edge of our fi eld of view. The 6 AU belt seen by 


( 2004|) is below our inne r working angle. 
'Ukamoto et al.l (|2004|) belt at ^16 AU 


tie 


Okamoto et al. 

We note that t 
only occurs on the INF side, at roughly the same location 
as the tentative brightness asymmetry seen in our disk 
model residuals. We see no evidence of the other belts in 
projection, but we model the disk with a continuous dust 
distribution and therefore may not b e sensitive to dust 
at the ir locations. Mid-IR imaging by [Weinberger et al 


(2003) indicates emission within 20 AU that is signih- 


cantly offset in position angle from the main outer disk. 
In our disk image we see no indication of this offset. 

Previous studies of (3 Pic’s debris disk in polarized 
scat tered light have been carried out both in the opti¬ 
cal (iGle dhill et al.jjlQ^ [Wo lstencroft et al.||1995|) and 
the NIR (jTamura et al.|2006p. These observations image 
the disk at separations of 15" — 30" and 2.6" — 6.4", in 
the optical and NIR, respectively. At these angular sep¬ 
arations the total intensity observations are not limited 
by the PSF halo and when combined with the polar¬ 
ized images, polarization fraction can be used to model 








































































Figure 9. Three different disk models with polarized intensity contours overplotted. The differences in the shapes and spacing of the 
contours illustrate how the inner and outer radius are constrained even though neither are directly detected. All three images are displayed 
with the same Log colour scale. This colour scale has been chosen to em^asize the differences between the models, and is not the same 
scale as Figure^ and Figure Left: The best fit disk model from Figure^ Center: The same disk model with the inner radius changed 
from the median value of Ri = 23.6 AU to 15 AU. The smaller inner radius increases the scattering contributions from dust at smaller 
projected separations. As a result, the spacing of the inner contours, in particular in the horizontal direction, becomes tighter while leaving 
the outer contours relatively unchanged. In addition, the contour ansae are pulled down towards midplane. Right: The best fit disk model 
with the outer radius changed from R 2 = 138 AU to 200 AU. The larger outer radius increases the scattering contributions from dust at 
separations further above the midplane, which pushes the contours further out in the vertical direction. 

within 0.05° from 90°. Th is is a reduction by a fac tor of 
^ 50 from the estimate in Macintosh et al. (2014), who 
found i = 90.7 ± 0.7. 


the dust g:r ai ns ([Voshchinnikov fc Krug:el||1999| |Krivova 
et al.||2000|). |T'amura et al.| (|2006p combine the optical 
measurements with their A-band data and find that the 
observations could be explained by scattering from fluffy 
aggregates made up of sub-micron dust grains. Unfor¬ 
tunately, the lack of total intensity images and a non¬ 
overlapping field of view make a direct comparison be¬ 
tween our observations and this past work difficult. 

5.2. P Pic b 


viously published (e.g. 

Chauvin et al. 

2012 

et al.||2014 Nielsen et a 

i.|2014), but the long( 


Macintosh 


tighten the constraints on the orbital parameters. In par¬ 
ticular, we find that the position angle of the ascending 
node of the planet lies in betwee n the main outer di sk 


and warp feature, consistent with Nielsen et al. (2014). 
At first glance, the errors o n our orbital elements ap 


pear comparable to those in Macintosh et al. (2014) 


However, our fit includes the total mass of the Ast ern as 


an additional free parameter. Nielsen et al. (2014) mod¬ 
eled the system’s total mass as a free parameter m their 
orbital fit and found that accounting for the uncertainty 
in the system’s total mass resulted in larger uncertainties 
in the planet’s orbital elements. In particular, they find 
that with a floating system mass the eccentricity distri¬ 
bution has a long tail that peaks at high eccentricities. 
Due to a degeneracy between semi-major axis and eccen¬ 
tricity, this stretched the semi-major axis distribution to 
higher values as well. In Figurewe And that the eccen¬ 
tricity is now significantly better constrained (e < 0.26), 
and while the degeneracy remains, the semi-major axis 
is constrained to be < 10.7 AU with 84% confidence. 

For each orbit defining our posterior distribution, we 
calculate the epoch of closest approach and And that it 
will fall between 2017 November 20, and 2018 April 4 
with 68% confidence. With our derived inclination of i = 
89.01±0.36°, the updated transit probability is ^ 0.06%, 
assuming that the planet will transit if the inclination is 


Even though the likelihood of a planet transit is small, 
it is still possible that dust particles orbiting within the 
planet’s Hill sphere {Rmw ~ 1 AU) will transit. Indeed 
the transit of a ring system surrou nding an exoplanet wa s 


recently detected around J1407 ([Mamajek et al. 2012). 


In the outer solar system, satellites around the giant 
planets have stable orbits within a Hill sphere about the 
planet out to ^ O.bRum whe n in prograde orbits and 
0.7i?Hiii in retrograde orbits ( Shen fc Tremain^|2008 ). 


For P Pic 6, assuming a planetary mass of 11 IVlj^ a 
semimajor axis of 9.25 AU, a circular orbit and a stellar 
mass of 1.61 Mq, we calculate a Hill radius of ^ 1.2 AU. 
Thus, stable orbits within the Hill sphere will transit if 
the planet’s inclination is within 3.8° and 5.3° of edge-on, 
for prograde and retrograde orbits, respectively. Our new 
constraints on the inclination indicate that these orbits 
will almost certainly transit. However, the true transit 
probability will depend not only on the exact location of 
the dust, but also its orientation relative to the observer. 
For example, dust that Alls the stable orbits and is or¬ 
biting face-on relative to the observer will transit, but if 
it is orbiting edge-on it will not. 

The presence of infalling comets (a.k.a. falling evapo¬ 
rating bodies, or FEBs) has been previously infer red by 


& Morbidelli| 

fi^ 

). Thebault & Beustj 

1 (2001) 

that a massive [M 

> iUjup) planet witJ 

im ^ 2 


slightly eccentric orbit (e > 6.05 — 1), could be responsi¬ 
ble for imparting highly elliptical orbits on bodies within 
a 3:1 or 4:1 resonance, that then plunge towards the star. 
In this scenario the argument of the periastron of the 
planet is restricted to a value oi uo' = — 70 ± 20° from 
the line of sight. Using our deflnitions, the equivalent 
requirement is cj = 200 ± 20°. Our results neither con- 
Arm nor rule out the infalling comet scenario. While the 
marginalized distribution of the argument of periastron 
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allows for a broad range of values, if the orbit is indeed 
eccentric, then uj peaks strongly around 170°, just out¬ 
side of the acce ptable values for this scenario. [Thebault] 
& Beust| (|2001|) find that if the eccentricity of the mas¬ 
sive perturber (here assumed to be (3 Pic h) is as large 
as e ~ 0.1 then the infalling comets most likely originate 
in the 3:1 resonance, which occurs between 18 AU and 
22 AU based on our 68% confidence range for P Pic h. 
The inner edge of the dust in our scattered light images 
falls at Ri = 23.6llo.57 AU, outside of range of values for 
the 3:1 resonance. However, as noted above, our inner 
radius is sharply defined, and there may still be material 
inside. For smaller perturber eccentricities the infalling 
comets originate in the 4:1 resonance, which occurs be¬ 
tween 22.2 AU and 26.5 AU. Our disk model does not 
constrain whether there is an excess of bodies librating 
in the 4:1 resonance. 

5.3. The disk-planet interaetion 

A planet on an inclined orbit is thought to be re¬ 
sponsible for the warp feature in the re gion of the disk 
outside our field of view a t ^ 80 AU (Mouillet et al. 
1997 An gereau et al.||2QQl|). T he directly-imaged planet 
P Pic b ([Lagrange et al. 2009) has a mass, semi-major 


axis, and inclination cons istent with producing the warp 


(e.g., Dawson et al. ]MTt . The updated position angle of 
P Pic 6's ascending node, U = 31.75 ± 0.15°, is offs et by 
2.65° with respect to the flat outer disk (29.1±0.1°;|Apai 
et al. 2015), consistent with producing a warp tilted by 5"" 
counter-clock wise with respect t o the flat outer disk. As 
illustrated by |Apai et al.|[?015[ Fig. 21, our azimuthal 
viewing angle of the warped disk affects the degree to 
which the inner disk and the planet’s orbit appear aligned 
with the flat outer disk in projection and also affects the 
projected height of the warp. Although the planet’s up¬ 
dated orbit remains consistent with sculpting the outer 
regions of disk, several features of the inner regions of 
the disk that we measured here are unexpected solely 
from sculpting by P Pic b (Figure [Tq|). 

First, the inner edge of the disk is at 23.6 AU, about 
twelve Hill radii from t he planet. We pe rformed a simula¬ 
tion using mercury6 ( Chambers! 1999) of a planet with 
orbital parameters set to the median values in Figure 3 
embedded in a disk of test particles initially spanning 10 
to 40 AU. On the timescale of hundreds to thousands of 
orbits, the planet clears out the disk to ^ 15 AU, with 
the inner edge p ersisting at that location over the 20 Myr 


stella r lifetime (Binks & Jeffries 2014 Mamajek & Bell 
2014). An inn er edge at ^ 15 AU i s m agreement with 
simulations by Rodigas et al. (2014|); cf. their Table 2. 

We have not explored disk models with gradual inner 
edges (e.g., Milli et al. 2014), so there may be material 
between 15 and 23.6 AU with a lower surface density, 
or planet bodies that are less collisionally active. If the 
disk inside 23.6 AU is truly cleared out, an undetected 
low-mass planet in between P Pic b and the disk’s inner 
edge could be responsible; we find that a planet could 
exist on a stable orbit in that region. 

Second, we expect the inner disk to be centered on 
the planet’s orbital plane. Given a warp located at 
^ 80 AU, the width of a secular cycle (i.e., the differ¬ 
ence in semi-major axis for which the planetesimals are 
27r out of phase in their oscillations about the planet’s 
orbital plane) is only about 1 AU at a radius of 40 AU 


and the timescale of a secular cycle is about forty times 
shorter than at the location of the warp. Therefore, 
close to the planet, a sufficient number of secular cy¬ 
cles should have passed that the parent bodies’ free in¬ 
clination vectors are randomized about the forced in¬ 
clination from the planet. Under certain conditions, 
we found that our simulation could produce a parent 
bodies sky plane inclinations distribution with peaks at 
^ “^planet,sky ipianet,outerdisk (onc of which could Corre¬ 
spond to ^ 85°), where ipianet,sky is the line of sight in¬ 
clination of the planet and ipianet,outerdisk is the mutual 
inclination between the planet and the outer disk. How¬ 
ever, we expect that even in these circumstances the mea¬ 
sured disk midplane would be aligned with the planet’s 
orbital plane. Moreover, damping by collisions, small 
bodies, or residual gas—provided that it occurs on a 
timescale shorter than half a secular timescale—reduces 
the free inclination, decreasing the disk scale height but 
keeping it centered about the planet’s orbital plane. 

Instead, the average plane of the inner disk appears 
mutually inclined with respect to the planet’s orbit. If 
the polarized intensity images were dominated by scat¬ 
tered light from the outer disk, a mutual inclination with 
respect to the planet could be consistent, (depending on 
the s emi-major axis of th e dominant dust; see Figure 1 
from Dawson et al.||2011), but in the current disk model 
the observed light is dominated by a close-in disk. Con¬ 
tribution from another planet to the forced plane of the 
disk is a possibility but the available parameter space 
for an additional planet that tilts the disk toward us, 
yet is too low mass to escape detection, is quite lim¬ 
ited. In the future, we plan to explore a wider range of 
dust-scattering models to ensure that this result (a disk 
mutually inclined to the planet’s orbit at ^ 25 AU) is 
not dependent on the assumed dust properties. 

Finally, the scale height of the disk appears larger than 
expected from stirring by p Pic b or self-stirring of the 
parent bodies. In the absence of damping, the total 
thickness of the disk would be ^ 2ip, corresponding to 
a scale height aspect ratio of about 0.06 for a planet in¬ 
clined by 3.6° with respect to the primordial plane. Self¬ 
stirring to the escape velocity of 10 km planetesimals 
would contribute only about 0.001 to the aspect ratio; 
self-stirring to the escape velocity of Pluto sized bodies 
would be required. In practice, we do not expect most 
parent bodies participating in the collisional cascade to 
be stirred to random velo cities of the largest bodies (e.g.. 


Pan & Schlichting 2012); their steady-state random ve- 
locities depend on the balance between stirring, damping 
by smaller bodies and each other, and radiation forces. 
The scale height is also si gnificantly larger tha n observed 
further out in the disk (Ahmic et al. 2009)—even at 
50 AU dMilli et al.||2014|). 'I'he robustness of the scale 
height to the dust scattering model should be explored 
further; for example, a significant contribution from po¬ 
larized back scattering (not modeled here) could result 
in a smaller inferred scale height. If the current inferred 
large scale height in the very inner disk is robust, a sub¬ 
detection planet located between P Pic b and the inner 
edge of the disk and mutually inclined with respect to 
P Pic 6 is a possible explana tion. 


Nesvold & Kuchner (2015) recently simulated the dy- 
namical and collisional beha vior of p Pic’s planet esimals 
and dust grains using SMACK ( Nesvold et al.|2013 ), which 
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Figure 10. Three orthogonal projections of the system: the plane of the sky {bottom left), a top down view {top left), and a side view 
{bottom right). The system has been rotated in three dimensions so that the midplane of outer disk is horizontal (blue line) in the bottom 
left plot. Each image includes the best fit disk midplane (greyscale decreasing as r“^), a random selecti on of 100 accepted orb i ts (dotted red 
lines) and the location of /3 Pic b according to a likely orbit at the same epochs as the measurements of|Chauvin et al.|(|2012|),|Nielsen et al.| 
(|2014|) and those included in this work (purple markers). Bottom left: The positions of the planet emulate the direct imaging astrometry 
points. The green line indicates the position angle of our inner disk model. Top left: The orbital inclination of the planet and the mismatch 
between the planet and main disk’s position angles on the sky result in the top panel being slightly offset from a planet’s orbital plane. 
The RV measurement of the planet ([Snellen et al. |2014t breaks the degeneracy in the orbital direction and allows us to calculate the line 
of sight {Z) coordinate for each epocn. i he red arrow indicates the direction of motion in the orbit. Bottom right: The green line displays 
the inclination of the inner disk relative to the observer’s point of view. In all panels, the green x indicates the location of the star. 
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models planetesimals across a range of sizes using super 
particles. They find that collisional damping is not im¬ 
portant in shaping the morphology of the disk. Their 
detailed model also does not predict the surprising ob¬ 
servational features discovered here: they find the planet 
clears a gap only out to 14.5 AU and that the disk is cen¬ 
tered about the planet’s orbital plane (see their Figure 
3). They find that some planetesimals in the inner disk 
are scattered by each other or the planet to inclinations 
larger than increasing the thickness of the inner disk 
by about 50%, not enough to account for the 200% 
larger) observed scale height. 

6. CONCLUSION 

We have presented new images of the P Pic debris disk 
in polarized light that reach angular separations previ¬ 
ously inaccessible to both space and ground-based tele¬ 
scopes. The use of PDI as a means of PSF subtraction 
circumvents the need for ADI PSF subtraction which can 
cause self-subtraction, especially in vertically extended 
disks like that of /3 Pic at the angular separations ex¬ 
plored by GPL The disk image was modeled with a radial 
power-law dust distribution combined with a Henyey- 
Greenstein scattering function. The disk model indicates 
an inclined disk at a position angle on the sky between 
the main outer disk and the warped feature with an inner 
edge at ^ 23 AU. 

The conclusions about the geometry of the disk are 
based on the assumption that a Henyey-Greenstein scat¬ 
tering phase function can accurately represent the true 
scattering properties of the constituent dust grains. Fu¬ 
ture imaging studies, such as multi-color polarimetry at 
similar angular separations, will allow for the use of more 
sophisticated dust grain models that will be able to fur¬ 
ther examine the inner part of the disk and to test our 
results. 

In addition, we presented ten new astrometric mea¬ 
surements of the planet (3 Pic 6, which we combine with 
previous measurements to fit an orbital solution. The 
solution improves upon those previously published by 
tightening the constraints on the Keplarian orbital el¬ 
ements, particularly the inclination and position angle 
of the ascending node. We leave the total mass of the 
system as a free parameter, allowing us to constrain the 
stellar mass of (3 Pic to within 5%. 

When considered together, the disk model and the or¬ 
bital fit indicate that the dynamics of the inner edge of 
the disk are not consistent with sculpting by the planet 
/d Pic 6 alone. This could be explained by an as-of-yet 
undetected planet in-between the known planet and the 
inner edge of the disk. Under this scenario the less mas¬ 
sive, further out planet would dynamically influence the 
inner regions of disk, while the more massive p Pic b 
would have a greater affect at larger radii, causing the 
well know warp. If there is in fact another planet at 
this location, this will have significant consequences for 
our understanding of the planet formation history and 
dynamical evolution of this system. 
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